#include "Signal.h"
#include <cassert>
using namespace std;

/*
// ---------------------------------------------------------------------
// CLASS: signal
// ---------------------------------------------------------------------
signal::signal(const char* iname,double it_end,double idt,int idimension,int ilog_mode)
{
int i,j;
string columnname;
char* filename;

name[0] = '\n';
strcpy(name,iname);
t_end = it_end;
dt = idt;
dimension = idimension;
log_mode = ilog_mode;

repeat = 0;
t_node = int(t_end/dt)+1;
time  = new double[t_node];
values = new double[dimension*t_node];
for(i = 0;i<t_node;i++)
{
time[i] = i*dt;
for(j=0;j<dimension;j++)
{
values[i*dimension+j] = 0.0;
}
}
t_idx = 0;

if(log_mode > 0)
{
means = new double[dimension*t_node];
vars  = new double[dimension*t_node];
for(i = 0;i<t_node;i++)
{
time[i] = i*dt;
for(j=0;j<dimension;j++)
{
means[i*dimension+j] = 0.0;
vars[i*dimension+j]  = 0.0;
}
}
}

if(log_mode == 2)
{
sprintf(filename,"%s_detail.dat",name);
log_detail.open(filename);
log_detail << "Signal: " << name;
log_detail << setw(COLUMNWIDTH)<< "\nTime ";
for(j=0;j<dimension;j++)
{
log_detail << setw(COLUMNWIDTH) << "signal[" << j << "]";
}
}
}
// ----------------------------------------------------------------------
signal::~signal()
{
int i,j;
char* filename;
char* columnname;

if(log_mode == 2)
{
log_detail.close();
}

if(log_mode>0)
{
filename = strcat(name,"_average.dat");
log_average.open(filename);
log_average << "Signal: " << name;
log_average << "\nRepeat = " << repeat;
log_average << setw(COLUMNWIDTH)<< "\nTime";
for(j=0;j<dimension;j++)
{
sprintf(columnname,"mean[%i]",j);
log_average << setw(COLUMNWIDTH) << columnname;
sprintf(columnname,"var[%i]",j);
log_average << setw(COLUMNWIDTH) << columnname;
}
for(i=0;i<t_node;i++)
{
log_average << "\n";
log_average << setw(COLUMNWIDTH) << time[i];
for(j=0;j<dimension;j++)
{
log_average << setw(COLUMNWIDTH) << means[i*dimension+j];
log_average << setw(COLUMNWIDTH) << vars[i*dimension+j];
}
}
log_average.close();
delete []means;
delete []vars;
}

delete []time;
delete []values;
}
// ----------------------------------------------------------------------
void signal::reset()
{
repeat++;
t_idx = 0;
}
// ----------------------------------------------------------------------
void signal::set(double ivalue)
{
value[t_idx] = ivalue;
t_idx++;
}
// ----------------------------------------------------------------------
void signal::set(double *ivalue)
{
int j;
for(j=0;j<dimension;j++)
{
value[t_idx*dimension+j]=ivalue[j];
}
}
// ----------------------------------------------------------------------
double signal::get(double itime)
{
return values[t_idx*dimension];
}
// ----------------------------------------------------------------------
*/




// ----------------------------------------------------------------------
//
// ClASS: block
//
// ----------------------------------------------------------------------
block::block()
{
#ifdef DEBUG
  cout << "block::block()\n";
#endif
  time = 0.0;

  input_num = 0;
  input_dim = NULL;
  inputs = NULL;

  output_num = 0;
  output_dim = NULL;
  outputs = NULL;
  // outputs_next = NULL;
#ifdef DEBUG
  cout << "END block::block()\n";
#endif
}
// ----------------------------------------------------------------------
//block::block(const block & blk){
//  cout << "block::block(const block & blk)\n";
//  int i,j;
//  time = 0.0;
//
//  input_num = 0;
//  input_dim = NULL;
//  inputs = NULL;
//
//  output_num = 0;
//  output_dim = NULL;
//  outputs = NULL;
//  outputs_next = NULL;
//
//  initialize(blk.input_num,blk.input_dim,
//    blk.output_num,blk.output_dim,blk.dt);
//  time = blk.time;
//  for(i=0;i<input_num;i++){
//    for(j=0;j<input_dim[i];j++){
//      inputs[i][j] = blk.inputs[i][j];
//    }
//  }
//  for(i=0;i<output_num;i++){
//    for(j=0;j<output_dim[i];j++){
//      outputs[i][j] = blk.outputs[i][j];
//      outputs_next[i][j] = blk.outputs_next[i][j];
//    }
//  }
//  cout << "END block::block(const block & blk)\n";
//}
//// ---------------------------------------------------------------------
//block & block::operator=(const block & blk){
//  cout << "block::operator=()\n";
//  if(this==&blk){
//    return *this;
//  }
//  int i,j;
//  if(input_num>0){
//    delete [] input_dim;
//    delete [] inputs;
//  }
//  if(output_num>0){
//    delete [] output_dim;
//    for(i=0;i<output_num;i++){
//      delete [] outputs[i];
//      delete [] outputs_next[i];
//    }
//    delete [] outputs;
//    delete [] outputs_next;
//  }
//  initialize(blk.input_num,blk.input_dim,
//    blk.output_num,blk.output_dim,blk.dt);
//  time = blk.time;
//  for(i=0;i<input_num;i++){
//    for(j=0;j<input_dim[i];j++){
//      inputs[i][j] = blk.inputs[i][j];
//    }
//  }
//  for(i=0;i<output_num;i++){
//    for(j=0;j<output_dim[i];j++){
//      outputs[i][j] = blk.outputs[i][j];
////      outputs_next[i][j] = blk.outputs_next[i][j];
//    }
//  }
//  cout << "END block::operator=()\n";
//  return *this;
//}
// ---------------------------------------------------------------------
block::block(int i_num,const int* i_dim,
	     int o_num,const int* o_dim,
	     double idt)
{
  cout << "block::block(int i_num, const int* i_dim, int o_num, const int* o_dim, double idt)\n";
  initialize(i_num,i_dim,o_num,o_dim,idt);
  cout << "END block::block(int i_num, const int* i_dim, int o_num, const int* o_dim, double idt)\n";
}
// ----------------------------------------------------------------------
block::~block()
{
#ifdef DEBUG
  cout << "block::~block()\n";
#endif

  int i;
  if(input_num>0){
    delete []input_dim;
    delete []inputs;
  }
  if(output_num>0){
    delete []output_dim;
    for(i=0;i<output_num;i++)
    {
      delete [] outputs[i];
//      delete [] outputs_next[i];
    }
    delete []outputs;
//    delete []outputs_next;
  }

#ifdef DEBUG
  cout << "END block::~block()\n";
#endif
}
// ----------------------------------------------------------------------
bool block::initialize(int i_num, const int *i_dim, 
		       int o_num, const int *o_dim, 
		       double idt)
{
#ifdef DEBUG
  cout << "block::initialize()\n";
#endif
  
  int i,j;
  input_num = i_num;
  if(input_num>0){
    input_dim = new int[input_num];
    if(input_dim==NULL){
      cout << "New failed\n";
      exit(1);
    }
    for(i=0;i<input_num;i++){
      input_dim[i] = i_dim[i];
    }
    inputs = new double*[input_num];
    for(i=0;i<input_num;i++){
      inputs[i] = NULL;
    }
  }
  // The block is a source
  else{
    input_dim = NULL;
    inputs = NULL;
  }

  output_num = o_num;
  if(output_num>0){
    output_dim = new int[output_num];
    for(i=0;i<output_num;i++){
      output_dim[i] = o_dim[i];
    }
    outputs = new double*[output_num];
    //outputs_next = new double*[output_num];
    for(i=0;i<output_num;i++){
      outputs[i] = new double[output_dim[i]];
      //outputs_next[i] = new double[output_dim[i]];
      for(j=0;j<output_dim[i];j++){
	outputs[i][j] = 0.0;
	//outputs_next[i][j] = 0.0;
      }
    }
  }
  // The block is a sink
  else{
    output_dim = NULL;
    outputs = NULL;
    //outputs_next = NULL;
  }

  time = 0;// NOTE: For dynamic system, the starting time is dt rather than 0.
  dt = idt;
#ifdef DEBUG
    cout << "input_num=" << input_num << "\n";
    cout << "input_dim=" << input_dim << "\n";
    for(i=0;i<input_num;i++){
      input_dim[i] = i_dim[i];
      // cout << "input_dim = " << input_dim << "\n";
      cout << "input_dim[" << i << "]=" << input_dim[i] << "\n";
    }
    cout << "END block::initialize()\n";
#endif 
  return true;
}
// ----------------------------------------------------------------------
bool block::ConnectInputTo(int input_ID,block& source,int output_ID)
{
  assert(input_ID<input_num);
  assert(output_ID<source.getOutputNum());
  // assert(input_dim[input_ID]==source.getOutputDim(output_ID));

  if(input_dim[input_ID]!=source.getOutputDim(output_ID)){
    cerr << "input_dim" << input_dim << "\n";
    cerr << "input_dim[" << input_ID << "]=" << input_dim[input_ID] << "\n";
    cerr << "source.getOutputDim(" << output_ID << ")="<< source.getOutputDim(output_ID) << "\n";
    exit(0);
  }
  inputs[input_ID] = source.getOutputAddress(output_ID);
  return true;
}
//// ----------------------------------------------------------------------
//void block::reset()
//{
//  time = 0;// NOTE: For dynamic system, the starting time is dt rather than 0.
//}
//// ----------------------------------------------------------------------
//void block::output()
//{
//  int i,j;
//  for(i=0;i<output_num;i++){
//    for(j=0;j<output_dim[i];j++){
//      outputs[i][j] = outputs_next[i][j];
//    }
//  }
//}
/*
// ----------------------------------------------------------------------
TimesTwo::TimesTwo(double idt):block("TimesTwo",1,1,idt)
{
}
// ----------------------------------------------------------------------
void TimesTwo::update(double itime,double idt)
{	
assert(itime >= time[t_idx]);

if(itime+idt<time<time[t_idx]+dt)
{
return; // Do nothing
}
else{itime+idt>time+dt) // Update the output
{
I1 = inputs[1]-->getValue();
O1 = I1*2;
outputs[1]-->setValue(O1);
}
}
*/
// ----------------------------------------------------------------------



// ----------------------------------------------------------------------
//
// class: simulation
// 
// ----------------------------------------------------------------------
simulation::simulation(double idt, double it_end, int iRepeat, int iJobID, int iMode){
  dt     = idt;
  t_end  = it_end;
  repeat = iRepeat;
  JobID  = iJobID;
  Mode   = iMode;
  time(&time_begin);
}

simulation::~simulation(){	
  time(&time_end);
  ofstream flog("simulation.log");
  flog << "INPUTS ===============================\n";
  flog << "Job ID  : " << JobID << "\n";
  flog << "t_end(s): " << t_end << "\n";
  flog << "dt      : " << dt << "\n";
  flog << "Repeats : " << repeat << "\n";
  flog << "SUMMARY ==============================\n";
  flog << "Starts at: " << asctime(localtime(&time_begin));
  flog << "Ends at:   " << asctime(localtime(&time_end));
  flog << "Time(s):   " << difftime(time_end,time_begin);
  flog.close();
}
// ----------------------------------------------------------------------



// ----------------------------------------------------------------------
//
// class: MCSimulator
//
// ----------------------------------------------------------------------
MCSimulator::MCSimulator(int iEventNum){
  EventNum = iEventNum;
  Events = new int[EventNum];
  Prob   = new double[EventNum];

  int i;	
  for(i=0;i<EventNum;i++){
    Prob[i] = 0.0;
  }
  TotalProb = 0.0;

  // Copied from function time_seed() ----------
  time_t timeval;       /* Current time. */
  unsigned char *ptr;   /* Type punned pointed into timeval. */
  unsigned seed;        /* Generated seed. */
  size_t ii;

  timeval = std::time (NULL);
  ptr = (unsigned char *) &timeval;

  seed = 0;
  for (ii = 0; ii < sizeof timeval; ii++)
    seed = seed * (0xff + 2u) + ptr[ii];
  // -------------------------------------------

  gsl_rng_env_setup();
  rng = gsl_rng_alloc(gsl_rng_ranlxd1);
  gsl_rng_set(rng,seed);
}

MCSimulator::~MCSimulator(){
  delete [] Events;
  delete [] Prob;
}

int MCSimulator::TryOnce(){
  double x = gsl_rng_uniform(rng)*TotalProb;
  int i;
  double CummulatedProb = 0.0;
  for(i=0;i<EventNum;i++){
    CummulatedProb += Prob[i];
    if(CummulatedProb > x){break;} // Always use [a,b)(left close and right open) interval
  }
  return i;
}
/************************************************************************/




// ------------------------------------------------------------------
//
// CLASS: MCArray
//
// ------------------------------------------------------------------
MCArray::MCArray(int isize, gsl_rng *random_num_generator){
  int i,k;
#ifdef DEBUG
  cout << "MCArray::MCArray(" << isize << ")\n";
#endif
  size=isize;
  k = 0;
  while(pow(2.0,double(k))<size){
    k++;
  }// 2^(k)>=size
  array_size = int(pow(2.0,double(k+1))-1);
  array = new double[array_size];
  for(i=0;i<array_size;i++){
    array[i] = 0.0;
  }
  layer = k;
  array_index_begin = int(pow(2.0,double(k))-1);
  array_index_end   = array_index_begin+size-1;

  if(random_num_generator == NULL){
    EXTERNAL_RNG = false;
    gsl_rng_env_setup();
    gsl_T = gsl_rng_default;
    gsl_r = gsl_rng_alloc(gsl_T);
    gsl_rng_set(gsl_r,time_seed());    
  }
  else{
    EXTERNAL_RNG = true;
    gsl_r = random_num_generator;
  }

  // For debugging purpose
#ifdef DEBUG
  cout << "size= " << size << ", array_size=" << array_size << "\n";
  cout << "layer= " << layer << "\n";
  cout << "array_index_begin=" << array_index_begin << ", array_index_end=" << array_index_end << "\n";
  cout << "END MCArray::MCArray(int isize)\n";
#endif
}
// ------------------------------------------------------------------
MCArray::MCArray(const MCArray & ma){
#ifdef DEBUG
  cout << "MCArray::MCArray(const MCArray & ma)\n";
#endif
  size = ma.size;
  layer = ma.layer;
  array_size = ma.array_size;
  array_index_begin = ma.array_index_begin;
  array_index_end = ma.array_index_end;
  array = new double[array_size];
  int i;
  for(i=0;i<array_size;i++){
    array[i]=ma.array[i];
  }

  if(!EXTERNAL_RNG){
    gsl_rng_free(gsl_r);
  }

  cout << "END MCArray::MCArray(const MCArray & ma)\n";
}
// ------------------------------------------------------------------
MCArray & MCArray::operator=(const MCArray & ma){
  // cout << "MCArray & MCArray::operator=()\n";
  if(this==&ma)
    return *this;
  size = ma.size;
  layer = ma.layer;
  array_size = ma.array_size;
  array_index_begin = ma.array_index_begin;
  array_index_end = ma.array_index_end;
  delete []array;
  array = new double[array_size];
  int i;
  for(i=0;i<array_size;i++){
    array[i] = ma.array[i];
  }
  // cout << "END MCArray & MCArray::operator=()\n";
  return *this;
}
// ------------------------------------------------------------------
MCArray::~MCArray(){
  // cout << "MCArray::~MCArray()\n";
  // cout << "array address " << array << "\n";
  //  delete []array;
  delete []array;
  // cout << "array address after delete " << array << "\n";
  // cout << "END MCArray::~MCArray()\n";
}
// ------------------------------------------------------------------
void MCArray::SetElement(int index, double value){
  //  cout << "MCArray::SetElement(" << index << "," << value << ")\n";
  assert(index<size && index>=0);
  //	cout << "MCArray::SetElement()\n";
  // cout <<"array_index_begin=" << MCArray::array_index_begin << ",array_index_end=" << MCArray::array_index_end << "\n";
  // cout << "array[array_index_begin+index]=" << array[array_index_begin+index] << "\n";
  double delta = value-array[array_index_begin+index];
  if(delta==0.0){
    return;
  }
  // Update the array value
  array[array_index_begin+index]=value;
  // Update the parent nodes
  int i = array_index_begin+index;
  int k = layer;
  while(k>0){
    // i = int(floor((i-1)/2));
    i = parent(i);
    k--;
    // cout << "i="<<i<<", k="<<k<<"\n";
    // array[i] += delta;
    array[i] = array[left_child(i)]+array[right_child(i)];
  }		
  //	cout << "END of MCArray::SetElement()\n";
}
// ------------------------------------------------------------------
double MCArray::GetElement(int index){
  assert(index<size);
  return array[array_index_begin+index];
}
// ------------------------------------------------------------------
int MCArray::FindElement(double value){
  double value_tmp = value;
  int current_layer = 0;
  int i = 0;
  //for(i=0;i<int(pow(2.0,double(layer+1))-1);i++){
  //cout << "array["<<i<<"]="<<array[i]<<"\n";
  //}
  i = 0;
  while(current_layer<layer){
    if(value<array[left_child(i)]){
      i = left_child(i);			
    }	
    else{
      value = value-array[left_child(i)];
      i = right_child(i);			
    }
    current_layer++;
  }
  return i-array_index_begin;
}
// ------------------------------------------------------------------
int MCArray::PickAnEvent(){
  if(array[0]==0.0){
    for(int i=0;i<array_size;i++){
      cerr << array[i] << " ";
    }
    throw runtime_error("MCArray::PickAnEvent(): total probablity is zero");
  }
  double rdm = gsl_rng_uniform(gsl_r)*array[0];
  
  return FindElement(rdm);
}
// ------------------------------------------------------------------



/* ==================================================================
class: MathFunction
===================================================================*/
MathFunction::MathFunction(int i_num_input,int i_num_para){
  num_input = i_num_input;
  num_para = i_num_para;
}

MathFunction::~MathFunction(){
}

bool MathFunction::parameterize(const double *para){
  if(num_para!=0){
    cerr << "error @ MathFunction::parameterize()\n";
    return false;
  }
  return true;
}

bool MathFunction::parameterize(int para_ID,double para_val){
  if(num_para!=0){
    cerr << "error @ MathFunction::parameterize()\n";
    return false;
  }
  return true;
}

bool MathFunction::parameterize(int para_ID,double *para_val){
  return true;
}

void MathFunction::plot(const double *x_min,const double *x_max,const double *dx,string FileName){
  ofstream f_log(FileName.c_str());
  double f;
  switch(num_input){
  case 1:
    double x1[1];
    for(x1[0]=x_min[0];x1[0]<x_max[0];x1[0]=x1[0]+dx[0]){
      f = value(x1);
      f_log << setw(15) << x1[0] << setw(15) << f << "\n";
    }
    break;
  case 2:
    double x2[2];
    for(x2[0]=x_min[0];x2[0]<x_max[0];x2[0]=x2[0]+dx[0]){
      for(x2[1]=x_min[1];x2[1]<x_max[1];x2[1] = x2[1]+dx[1]){
	f = value(x2);
	f_log << setw(15) << x2[0] << setw(15) << x2[1] << setw(15) << setprecision(10) << f << "\n";
      }
    }
  }
}

// ----------------------------------------------------------------
int perturb_MathFunction(MathFunction* fun,const double *x,const double *epsilon){
  cout << "***********************************************\n";
  cout << "Verfication of gradient calculation\n";
  int num_input = fun->get_num_input();
  double *grad_ana = new double[num_input];// analytical expression of gradient
  double *x_perturbed = new double[num_input];
  double val;
  int i,j;
  int RETURN_SIGN=0;
  cout << "x = [" << x[0];
  for(i=1;i<num_input;i++){
    cout << ", " << x[i];
  }
  cout << "]\n";
  fun->gradient(x,grad_ana);
  val=fun->value(x);
  cout << "val witout pert.:" << val << "\n";
  // Finite difference approximation

  double grad_app_pos,grad_app_neg;
  double val_pos,val_neg;
  for(i=0;i<num_input;i++){
    // Positive perturbation
    for(j=0;j<num_input;j++){
      if(j==i){
	x_perturbed[j] = x[j]+epsilon[j];
      }
      else{
	x_perturbed[j] = x[j];
      }
    }
    val_pos = fun->value(x_perturbed);
    grad_app_pos = (val_pos-val)/(epsilon[i]);


    // Negative perturbation
    for(j=0;j<num_input;j++){
      if(j==i){
	x_perturbed[j] = x[j]-epsilon[j];
      }
      else{
	x_perturbed[j] = x[j];
      }
    }
    val_neg = fun->value(x_perturbed);
    grad_app_neg = (val_neg-val)/(-epsilon[i]);


    cout << "Grad. w.r.t x[" << i << "], perturbation = +/-(" <<epsilon[i] << ")\n";
    cout << "analytical exp: " << grad_ana[i] << "\n";
    cout << "val with pos. pert.:" << val_pos << "\n";
    cout << "positive pert.: " << grad_app_pos <<"\n";
    cout << "val with neg. pert.:" << val_neg << "\n";
    cout << "negative pert.: " << grad_app_neg << "\n";

    //    if(!((grad_app_pos<=grad_ana[i] && grad_ana[i]<=grad_app_neg)||
    // (grad_app_neg<=grad_ana[i] && grad_ana[i]<=grad_app_pos))){
    double grad_app_cen = (grad_app_pos+grad_app_neg)/2;
    if(abs((grad_app_cen-grad_ana[i])/grad_ana[i])>0.05){
      RETURN_SIGN = -1;
    }
  }
  delete []grad_ana;
  delete []x_perturbed;
  return RETURN_SIGN;
}
// ----------------------------------------------------------------
void MathFunction_plot(MathFunction* fun, 
		      const double *x_min,
		      const double *x_max,
		      const double *dx,		      
		      const int dimension,
		      bool (*FuncPtrValidRegime)(double* x,int dimension),
		      const double VALUE_FOR_INVALID_REGIME,
		      string FileName){
  string FileName_x1 = FileName+"_x1.dat";
  string FileName_y = FileName+"_y.dat";
  ofstream log_x1(FileName_x1.c_str());
  ofstream log_y(FileName_y.c_str());
  double x[2];

  switch(dimension){
    case 1:
      for(x[0]=x_min[0];x[0]<x_max[0];x[0]+=dx[0]){
	log_x1 << setw(15) << x[0] << "\n";
	log_y  << setw(15) << fun->value(x) << "\n";
      }    
      break;
    case 2:
      string FileName_x2 = FileName+"_x2.dat";
      ofstream log_x2(FileName_x2.c_str());
      for(x[0]=x_min[0];x[0]<x_max[0];x[0]+=dx[0]){
	log_x1 << setw(15) << x[0] << "\n";
      }
      for(x[1]=x_min[1];x[1]<x_max[0];x[1]+=dx[1]){
	log_x2 << setw(15) << x[1] << "\n";
      }
      // Output format: row --> x2, column --> x1
      // Example:
      //        x1[0]   x1[1]
      // x2[0]  y       y
      // x2[1]  y       y
      for(x[1]=x_min[1];x[1]<x_max[1];x[1]+=dx[1]){
	for(x[0]=x_min[0];x[0]<x_max[0];x[0]+=dx[0]){
	  if((*FuncPtrValidRegime)(x,2)){
	    log_y << setw(15) << fun->value(x);
	  }
	  else{
	    log_y << setw(15) << VALUE_FOR_INVALID_REGIME;
	  }
	}
	log_y << "\n";
      }
      log_x2.close();
      break;
  }
  log_x1.close();
  log_y.close();
}

/************************************************************************

RANDOM NUMBER GENERATION

************************************************************************/
/* Choose and return an initial random seed based on the current time.
Based on code by Lawrence Kirby <fred@genesis.demon.co.uk>.
Usage: srand (time_seed ()); */
unsigned time_seed (void)
{
  time_t timeval;	/* Current time. */
  unsigned char *ptr;	/* Type punned pointed into timeval. */
  unsigned seed;	/* Generated seed. */
  size_t i;

  timeval = time (NULL);
  ptr = (unsigned char *) &timeval;

  seed = 0;
  for (i = 0; i < sizeof timeval; i++)
    seed = seed * (0xff + 2u) + ptr[i];

  return seed;
}
/***********************************************************************/
double ranO(int *idum)
/*Returns a uniform random deviate between 0.0 and 1.0 using rand().*/
/*Set idum to any negative value to initialize or reinitialize the sequence*/
{
  static double y, maxran, v[131];
  double dum;
  static int iff=0;
  int j;

  if (*idum < 0 || iff == 0)
  {
    iff = 1;
    maxran = RAND_MAX+0.0;
    srand((unsigned)time(NULL));
    *idum = 1;
    for (j=0; j<=130; j++)
    {
      dum=rand()+0.0;
    }
    for (j=0; j<=130; j++)
    {
      v[j]=rand()+0.0;
    }
    y=rand()+0.0;
  }
  j= (int)(130.0*y/maxran);
  if (j > 130 || j < 0)
  {
    printf("RanO: This cannot happen.\n");
  }
  y=v[j];
  v[j]=rand();
  return y/(maxran+1.0);
}
/************************************************************************/
double RandRand()
{
  double r_tot, r;
  int i, n;
  static int seed;

  n=3;
  r_tot=0.0;

  for (i=0; i<n; i++)
  {
    r=ranO(&seed);
    if (i > 0)
    {
      r/=pow((RAND_MAX+0.0), (i+0.0));
    }
    r_tot+=r;
  }

  return r_tot;
}
/************************************************************************/
double RandN1()
{
  double rdm;

  while((rdm=RandRand())>=1.0) {}

  return rdm;
}
/************************************************************************/




/************************************************************************

OTHER FUNCTIONS

************************************************************************/
void xcov(double* x,double* xcov_val,int L){
  // Calculate auto-correlation function in the same way as Matlab function
  // XCOV
  // xcov(m) = \sum_{i=0}^{L-m-1}x[i]*x[i+m]
  int i,m;
  for(m=0;m<L-1;m++){
    xcov_val[m] = 0.0;
    for(i=0;i<L-1-m;i++){
      xcov_val[m]+=x[i]*x[i+m];
    }
  }
  return;
}
/************************************************************************/
